﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Threading;

using System.Numerics;

namespace Dwarf.Managed
{

    /// <summary>
    /// Solve the Spectral method problem.
    /// </summary>
    /// <remarks> 
    /// For the demonstration purposes of the project,
    /// we have chosen two-dimensional iterative FFT algorithm.
    /// See Solve() below
    /// </remarks>
    class Solver
    {
        /// <summary>
        /// Size of matrix.
        /// </summary>
        public int length;
        /// <summary>
        /// Working array.
        /// </summary>
        public Complex[][] complexArray;

        public Solver()
        {
        }

        /// <summary>
        /// Transposition of matrix
        /// </summary>
        private void transpose(Complex[][] matrix)
        {
            for (int i = 0; i < length; i++)
            {
                for (int j = i + 1; j < length; j++)
                {
                    var temp = matrix[i][j];
                    matrix[i][j] = matrix[j][i];
                    matrix[j][i] = temp;
                }
            }
        }

        private Complex[] W;
        /// <summary>
        /// 1-dimension Fast Fourier Algorithm
        /// </summary>
        private void FFT(Complex[] complexLine)
        {
            //variables for the FFT
            var length = complexLine.Length;
            int bitwidth = (int)Math.Log(length, 2);

            if (W == null || W.Length != length) {
                W = new Complex[length];

                //compute phase weights and reorder for butterflies
                for (var i = 0; i < length; i++) {
                    W[i] = new Complex(
                        Math.Cos(-((2 * i * Math.PI) / length)),
                        Math.Sin(-((2 * i * Math.PI) / length)));
                }
            }

            //reorder for butterflies
            int firstBit = length >> 1;
            int index = 0;
            int reverseIndex = 0;
            for (; ; ) {
                if (reverseIndex < index) {
                    var w = complexLine[index];
                    complexLine[index] = complexLine[reverseIndex];
                    complexLine[reverseIndex] = w;
                }
                index++;
                if (index == length)
                    // Break here, because the code below doesn't terminate for all 1's.
                    break;

                // Add 1 to reverse index, going from left to right:
                int carry = firstBit;
                // If adding 1 changed the bit to 0, we have a carry:
                while (((reverseIndex ^= carry) & carry) == 0)
                    carry >>= 1;
            }

            //combine and untangle applying phase weights
            //using successive butterflies of 2,4,...,bitwidth 
            for (var i = 1; i <= bitwidth; i++)
            {
                var halfStep = (1 << (i - 1));
                for (var j = 0; j < halfStep; j++)
                {
                    var tempShift = (j << (bitwidth - i));
                    var w = W[tempShift];
                    for (var k = j; k < length; k += halfStep * 2)
                    {
                        var u = complexLine[k + halfStep];
                        var v = complexLine[k];
                        var t = w * u;
                        complexLine[k + halfStep] = v - t;
                        complexLine[k] = v + t;
                    }
                }
            }
        }

        /// <summary>
        /// Based method.
        /// </summary>
        public void Solve()
        {
            /* Transform the rows */
            foreach (var row in complexArray) FFT(row);

            /* Transposition */
            transpose(complexArray);

            /* Transform the columns */
            foreach (var row in complexArray) FFT(row);

            /* Transposition */
            transpose(complexArray);
        }
    }
}
